#overwrites the zelig function so doesn't display legend
ci.plot <- function (obj, qi = "ev", var = NULL, ..., main = NULL, sub = NULL, 
          xlab = NULL, ylab = NULL, xlim = NULL, ylim = NULL, legcol = "gray20", 
          col = NULL, leg = 1, legpos = NULL, ci = c(80, 95, 99.9), 
          discont = NULL) 
{
  extract.sims <- function(obj, qi) {
    d <- length(obj$sim.out$range)
    k <- length(obj$sim.out$range[[1]][qi][[1]][[1]])
    hold <- matrix(NA, nrow = k, ncol = d)
    for (i in 1:d) {
      hold[, i] <- obj$sim.out$range[[i]][qi][[1]][[1]]
    }
    return(hold)
  }
  extract.sims1 <- function(obj, qi) {
    d <- length(obj$sim.out$range1)
    k <- length(obj$sim.out$range1[[1]][qi][[1]][[1]])
    hold <- matrix(NA, nrow = k, ncol = d)
    for (i in 1:d) {
      hold[, i] <- obj$sim.out$range1[[i]][qi][[1]][[1]]
    }
    return(hold)
  }
  ci.upper <- function(x, alpha) {
    pos <- max(round((1 - (alpha/100)) * length(x)), 1)
    return(sort(x)[pos])
  }
  ci.lower <- function(x, alpha) {
    pos <- max(round((alpha/100) * length(x)), 1)
    return(sort(x)[pos])
  }
  if (length(ci) < 3) {
    ci <- rep(ci, 3)
  }
  if (length(ci) > 3) {
    ci <- ci[1:3]
  }
  ci <- sort(ci)
  if ("timeseries" %in% obj$category) {
    xname <- "Time"
    qiseries <- c("pvseries.shock", "pvseries.innovation", 
                  "evseries.shock", "evseries.innovation")
    if (!qi %in% qiseries) {
      cat(paste("Error: For Timeseries models, argument qi must be one of ", 
                paste(qiseries, collapse = " or "), ".\n", sep = ""))
      return()
    }
    ev <- t(obj$getqi(qi = qi, xvalue = "x1"))
    d <- ncol(ev)
    xseq <- 1:d
    ev1 <- NULL
    if (is.null(xlab)) 
      xlab <- xname
    if (is.null(ylab)) {
      if (qi %in% c("pvseries.shock", "pvseries.innovation")) 
        ylab <- as.character(obj$setx.labels["pv"])
      if (qi %in% c("evseries.shock", "evseries.innovation")) 
        ylab <- as.character(obj$setx.labels["ev"])
    }
    if (is.null(main)) 
      main <- as.character(obj$setx.labels[qi])
    if (is.null(discont)) 
      discont <- 22.5
  }
  else {
    d <- length(obj$sim.out$range)
    if (d < 1) {
      return()
    }
    xmatrix <- matrix(NA, nrow = d, ncol = length(obj$setx.out$range[[1]]$mm[[1]]))
    for (i in 1:d) {
      xmatrix[i, ] <- as.matrix(obj$setx.out$range[[i]]$mm[[1]])
    }
    if (d == 1 && is.null(var)) {
      warning("Must specify the `var` parameter when plotting the confidence interval of an unvarying model. Plotting nothing.")
      return(invisible(FALSE))
    }
    xvarnames <- names(as.data.frame(obj$setx.out$range[[1]]$mm[[1]]))
    if (is.character(var)) {
      if (!(var %in% xvarnames)) {
        warning("Specified variable for confidence interval plot is not in estimated model.  Plotting nothing.")
        return(invisible(FALSE))
      }
    }
    if (is.null(var)) {
      each.var <- apply(xmatrix, 2, sd)
      flag <- each.var > 0
      min.var <- min(each.var[flag])
      var.seq <- 1:ncol(xmatrix)
      position <- var.seq[each.var == min.var]
    }
    else {
      if (is.numeric(var)) {
        position <- var
      }
      else if (is.character(var)) {
        position <- grep(var, xvarnames)
      }
    }
    position <- min(position)
    xseq <- xmatrix[, position]
    xname <- xvarnames[position]
    if (is.null(xlab)) 
      xlab <- paste("Range of", xname)
    ev1 <- NULL
    if (!is.null(obj$sim.out$range1)) {
      ev1 <- extract.sims1(obj, qi = qi)
    }
    ev <- extract.sims(obj, qi = qi)
    if (is.null(ylab)) {
      ylab <- as.character(obj$setx.labels[qi])
    }
  }
  k <- ncol(ev)
  n <- nrow(ev)
  if (is.null(col)) {
    myblue1 <- rgb(100, 149, 237, alpha = 50, maxColorValue = 255)
    myblue2 <- rgb(152, 245, 255, alpha = 50, maxColorValue = 255)
    myblue3 <- rgb(191, 239, 255, alpha = 70, maxColorValue = 255)
    myred1 <- rgb(237, 149, 100, alpha = 50, maxColorValue = 255)
    myred2 <- rgb(255, 245, 152, alpha = 50, maxColorValue = 255)
    myred3 <- rgb(255, 239, 191, alpha = 70, maxColorValue = 255)
    col <- c(myblue1, myblue2, myblue3, myred1, myred2, myred3)
  }
  else {
    if (length(col) < 6) {
      col <- rep(col, 6)[1:6]
    }
  }
  form.history <- function(k, xseq, results, ci = c(80, 95, 
                                                    99.9)) {
    history <- matrix(NA, nrow = k, ncol = 8)
    for (i in 1:k) {
      v <- c(xseq[i], median(results[, i]), ci.upper(results[, 
                                                             i], ci[1]), ci.lower(results[, i], ci[1]), ci.upper(results[, 
                                                                                                                         i], ci[2]), ci.lower(results[, i], ci[2]), ci.upper(results[, 
                                                                                                                                                                                     i], ci[3]), ci.lower(results[, i], ci[3]))
      history[i, ] <- v
    }
    if (k == 1) {
      left <- c(xseq[1] - 0.5, median(results[, 1]), ci.upper(results[, 
                                                                      1], ci[1]), ci.lower(results[, 1], ci[1]), ci.upper(results[, 
                                                                                                                                  1], ci[2]), ci.lower(results[, 1], ci[2]), ci.upper(results[, 
                                                                                                                                                                                              1], ci[3]), ci.lower(results[, 1], ci[3]))
      right <- c(xseq[1] + 0.5, median(results[, 1]), ci.upper(results[, 
                                                                       1], ci[1]), ci.lower(results[, 1], ci[1]), ci.upper(results[, 
                                                                                                                                   1], ci[2]), ci.lower(results[, 1], ci[2]), ci.upper(results[, 
                                                                                                                                                                                               1], ci[3]), ci.lower(results[, 1], ci[3]))
      v <- c(xseq[1], median(results[, 1]), ci.upper(results[, 
                                                             1], ci[1]), ci.lower(results[, 1], ci[1]), ci.upper(results[, 
                                                                                                                         1], ci[2]), ci.lower(results[, 1], ci[2]), ci.upper(results[, 
                                                                                                                                                                                     1], ci[3]), ci.lower(results[, 1], ci[3]))
      history <- rbind(left, v, right)
    }
    return(history)
  }
  history <- form.history(k, xseq, ev, ci)
  if (!is.null(ev1)) {
    history1 <- form.history(k, xseq, ev1, ci)
  }
  else {
    history1 <- NULL
  }
  if (k == 1) {
    k <- 3
  }
  all.xlim <- if (is.null(xlim)) 
    c(min(c(history[, 1], history1[, 1])), max(c(history[, 
                                                         1], history1[, 1])))
  else xlim
  all.ylim <- if (is.null(ylim)) 
    c(min(c(history[, -1], history1[, -1])), max(c(history[, 
                                                           -1], history1[, -1])))
  else ylim
  if (is.null(ylab)) 
    ylab <- "Expected Values: E(Y|X)"
  par(bty = "n")
  centralx <- history[, 1]
  centraly <- history[, 2]
  if (is.null(discont)) {
    gotok <- k
  }
  else {
    gotok <- sum(xseq < discont)
    if ((gotok < 2) | (gotok > (k - 2))) {
      cat("Warning: Discontinuity is located at edge or outside the range of x-axis.\n")
      gotok <- k
      discont <- NULL
    }
    if (gotok < k) {
      gotokp1 <- gotok + 1
      centralx <- c(centralx[1:gotok], NA, centralx[gotok + 
                                                      1:length(centralx)])
      centraly <- c(centraly[1:gotok], NA, centraly[gotok + 
                                                      1:length(centraly)])
    }
  }
  plot(x = centralx, y = centraly, type = "l", xlim = all.xlim, 
       ylim = all.ylim, main = main, sub = sub, xlab = xlab, 
       ylab = ylab)
  polygon(c(history[1:gotok, 1], history[gotok:1, 1]), c(history[1:gotok, 
                                                                 7], history[gotok:1, 8]), col = col[3], border = "white")
  polygon(c(history[1:gotok, 1], history[gotok:1, 1]), c(history[1:gotok, 
                                                                 5], history[gotok:1, 6]), col = col[2], border = "gray90")
  polygon(c(history[1:gotok, 1], history[gotok:1, 1]), c(history[1:gotok, 
                                                                 3], history[gotok:1, 4]), col = col[1], border = "gray60")
  polygon(c(history[1:gotok, 1], history[gotok:1, 1]), c(history[1:gotok, 
                                                                 7], history[gotok:1, 8]), col = NA, border = "white")
  if (!is.null(discont)) {
    polygon(c(history[gotokp1:k, 1], history[k:gotokp1, 1]), 
            c(history[gotokp1:k, 7], history[k:gotokp1, 8]), 
            col = col[3], border = "white")
    polygon(c(history[gotokp1:k, 1], history[k:gotokp1, 1]), 
            c(history[gotokp1:k, 5], history[k:gotokp1, 6]), 
            col = col[2], border = "gray90")
    polygon(c(history[gotokp1:k, 1], history[k:gotokp1, 1]), 
            c(history[gotokp1:k, 3], history[k:gotokp1, 4]), 
            col = col[1], border = "gray60")
    polygon(c(history[gotokp1:k, 1], history[k:gotokp1, 1]), 
            c(history[gotokp1:k, 7], history[k:gotokp1, 8]), 
            col = NA, border = "white")
    abline(v = discont, lty = 5, col = "grey85")
  }
  if (!is.null(ev1)) {
    lines(x = history1[1:gotok, 1], y = history1[1:gotok, 
                                                 2], type = "l")
    if (!is.null(discont)) {
      lines(x = history1[gotokp1:k, 1], y = history1[gotokp1:k, 
                                                     2], type = "l")
    }
    polygon(c(history1[1:gotok, 1], history1[gotok:1, 1]), 
            c(history1[1:gotok, 7], history1[gotok:1, 8]), col = col[6], 
            border = "white")
    polygon(c(history1[1:gotok, 1], history1[gotok:1, 1]), 
            c(history1[1:gotok, 5], history1[gotok:1, 6]), col = col[5], 
            border = "gray90")
    polygon(c(history1[1:gotok, 1], history1[gotok:1, 1]), 
            c(history1[1:gotok, 3], history1[gotok:1, 4]), col = col[4], 
            border = "gray60")
    polygon(c(history1[1:gotok, 1], history1[gotok:1, 1]), 
            c(history1[1:gotok, 7], history1[gotok:1, 8]), col = NA, 
            border = "white")
    if (!is.null(discont)) {
      polygon(c(history1[gotokp1:k, 1], history1[k:gotokp1, 
                                                 1]), c(history1[gotokp1:k, 7], history1[k:gotokp1, 
                                                                                         8]), col = col[6], border = "white")
      polygon(c(history1[gotokp1:k, 1], history1[k:gotokp1, 
                                                 1]), c(history1[gotokp1:k, 5], history1[k:gotokp1, 
                                                                                         6]), col = col[5], border = "gray90")
      polygon(c(history1[gotokp1:k, 1], history1[k:gotokp1, 
                                                 1]), c(history1[gotokp1:k, 3], history1[k:gotokp1, 
                                                                                         4]), col = col[4], border = "gray60")
      polygon(c(history1[gotokp1:k, 1], history1[k:gotokp1, 
                                                 1]), c(history1[gotokp1:k, 7], history1[k:gotokp1, 
                                                                                         8]), col = NA, border = "white")
    }
  }
  
  if (leg == 0)
    return()
  
  if (is.null(legpos)) {
    if (leg == 1) {
      legpos <- c(0.91, 0.04, 0.2, 0.05)
    }
    else if (leg == 2) {
      legpos <- c(0.09, 0.04, 0.2, 0.05)
    }
    else if (leg == 3) {
      legpos <- c(0.09, 0.04, 0.8, 0.05)
    }
    else {
      legpos <- c(0.91, 0.04, 0.8, 0.05)
    }
  }
  lx <- min(all.xlim) + legpos[1] * (max(all.xlim) - min(all.xlim))
  hx <- min(all.xlim) + (legpos[1] + legpos[2]) * (max(all.xlim) - 
                                                     min(all.xlim))
  deltax <- (hx - lx) * 0.1
  my <- min(all.ylim) + legpos[3] * min(max(all.ylim) - min(all.ylim))
  dy <- legpos[4] * (max(all.ylim) - min(all.ylim))
  lines(c(hx + deltax, hx + 2 * deltax, hx + 2 * deltax, hx + 
            deltax), c(my + 3 * dy, my + 3 * dy, my - 3 * dy, my - 
                         3 * dy), col = legcol)
  lines(c(hx + 3 * deltax, hx + 4 * deltax, hx + 4 * deltax, 
          hx + 3 * deltax), c(my + 1 * dy, my + 1 * dy, my - 1 * 
                                dy, my - 1 * dy), col = legcol)
  lines(c(lx - deltax, lx - 2 * deltax, lx - 2 * deltax, lx - 
            deltax), c(my + 2 * dy, my + 2 * dy, my - 2 * dy, my - 
                         2 * dy), col = legcol)
  lines(c(lx - 5 * deltax, lx), c(my, my), col = "white", lwd = 3)
  lines(c(lx - 5 * deltax, lx), c(my, my), col = legcol)
  lines(c(lx, hx), c(my, my))
  polygon(c(lx, lx, hx, hx), c(my - 3 * dy, my + 3 * dy, my + 
                                 3 * dy, my - 3 * dy), col = col[3], border = "white")
  polygon(c(lx, lx, hx, hx), c(my - 2 * dy, my + 2 * dy, my + 
                                 2 * dy, my - 2 * dy), col = col[2], border = "gray90")
  polygon(c(lx, lx, hx, hx), c(my - 1 * dy, my + 1 * dy, my + 
                                 1 * dy, my - 1 * dy), col = col[1], border = "gray60")
  polygon(c(lx, lx, hx, hx), c(my - 3 * dy, my + 3 * dy, my + 
                                 3 * dy, my - 3 * dy), col = NA, border = "white")
  text(lx, my, labels = "median", pos = 2, cex = 0.5, col = legcol)
  text(lx, my + 2 * dy, labels = paste("ci", ci[2], sep = ""), 
       pos = 2, cex = 0.5, col = legcol)
  text(hx, my + 1 * dy, labels = paste("ci", ci[1], sep = ""), 
       pos = 4, cex = 0.5, col = legcol)
  text(hx, my + 3 * dy, labels = paste("ci", ci[3], sep = ""), 
       pos = 4, cex = 0.5, col = legcol)
}